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Abstract 

We describe a novel method for obtaining the gradient expansion 
for the free energy of a clean BCS superconductor. We present explicit 
results up to fourth order in the gradients of the order parameter. 
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1 Introduction 



The Landau-Ginzburg equations [0 encapsulate much of the physics of con- 
ventional superconductors. They were originally proposed on phenomenolog- 
ical grounds, but after the advent of the BCS theory Gor'kov @ provided a 
formal derivation of the Landau-Ginzburg free energy, valid in the vicinity of 
T c . Gor'kov's derivation was soon extended to a wider range of temperature 
by Werthamer ||, and then to higher order in the gradients of the order 
parameter by Tewordt M and Eilenberger 

In the approach originally due to Gor'kov and then followed by Eilen- 
berger, Werthamer, Tewordt, and others, the free energy is calculated by 
relating its variation to the diagonal element of the Green function intro- 
duced by Gor'kov |§. These authors calculated the full Green's function in 
the Born series, and set its arguments equal at the end. Then they made an 
ansatz for the free energy in the same approximation by considering all possi- 
ble terms that may enter with undetermined coefficients. Upon the variation 
of the ansatz and comparison with the variation of the free energy, they ob- 
tained the coefficients. This is very laborious, and moreover nonsystematic. 

Since the super conducing gap is much smaller than the Fermi energy, it 
is usually safe to linearize the single-particle energy of the normal metal in 
the vicinity of the Fermi surface. Once this approximation is made, the full 
Gor'kov Green functions reduce to sums of Green functions, or resolvents, 
of one-dimensional differential operators 0. In the years since the publica- 
tion of the work described above, a great deal has been learned about the 
properties of such resolvents at points where the arguments coincide. When 
the Green function concerned is that of a Schrodinger operator, the diago- 
nal element of the resolvent satisfies an equation discovered by Gelfand and 
Dikii ||. This allows a fast and efficient evaluation of the diagonal element 
as an asymptotic expansion in powers of the potential and its derivatives. 
The terms in this expansion turn out to be the infinite sequence of con- 
served hamiltonian densities for the Korteweg-de Vries (KdV) hierarchy of 
integrable partial differential equations |J. 

For our superconductivity problem, the one-dimensional Schrodinger op- 
erator is replaced by a first-order matrix differential operator of the form 
studied by Andreev [p70| | . It is identical to a one-dimensional Dirac operator 
with a chiral mass term. Once again the diagonal element of the resolvent 
(or a simple modification of this) satisfies an equation which quickly and 
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efficiently generates a series expansion in the coeficients of the equation and 
their derivatives. The terms of this series give the hamiltonian densities 
for the family of integrable equations known as the AKNS hierarchy [|TT|], 
whose simplest member is equivalent to the non-linear Schrodinger equation 



Eq| . Perhaps surprisingly, this generating equation has also been exploited 
in the superconductivity literature. It is nothing other than the Eilenberger 
equation [13| which plays a central role in the quasiclassical theory [14 . 



In a previous paper JT5| we used both the Gelfand-Dikii equation and the 
Eilenberger equation to generate the gradient expansion for the free energy 
as far as the eighth order in the gradients of the gap function. Unfortunately 
these results are expressed in terms of gradients of the magnitude of the 
gap-function and gradients of its phase separately. This leads to very long 
expressions which are of limited utility. We were, however, able to expand the 
much more compact expressions appearing in the classic paper by Tewordt 
Pj] into our form and compare them. We found that our results differed from 
those of Tewordt at the fourth order — the highest order he had computed. 
While we had confidence in our own results (we had derived them from two 
completely different methods, using Mathematica to automate the labor), we 
were not able to isolate the source of the discrepancy. 

In this paper we report another method for deriving the gradient expan- 
sion. The new method is not quite as efficient at generating the series as 



that in ||15|| , but has the advantage that the resulting expressions are very 
compact. We evaluate all terms up to fourth order in the gradients, and 
expand them out so as to express our results in the form used by Tewordt 
[§]. We are thus able identify what appears to be a typographical error in 
his paper. While the correct fourth-order free energy is the principal result 
of the present paper, we believe that the methods we use are interesting in 
their own right, and will have application whenever one requires corrections 
to the leading orders of the quasiclassical method. 

In section two we briefly review how the calculation of the free energy of 
a three dimensional superconductor is reduced to computing determinants 
of one-dimensional differential operators. In section three we show how the 
Eilenberger equation arises as a property of the resolvent of such differential 
operators. In section four we use the Eilenberger function to derive a useful 
identity linking the dressing function for the Andreev hamiltonian to the 
determinant we wish to compute. In section five we show how this useful 
identity is related to the "shooting method" for computing determinants. 
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Then, in section six, we use our identity to compute the gradient expansion. 
Finally in section seven we compare our results with those of reference Q. 



2 From Three Dimensions to One 

In this section we will review standard material so as to establish our nota- 
tion. 

Following Gor'kov || we treat the normal metal as a free, highly degen- 
erate, electron gas and the superconductivity as arising from the BCS model 
potential - an instantaneous short-range attraction between pairs of elec- 
trons whose energy is within a narrow shell of width Wdebye about the Fermi 
surface. We also ignore paramagnetic effects. The partition function of such 
an electron gas may be written as a Berezin path integral 



Z = Tr (e- 0H ) 

d[ip]d[ip'] exp — 



d 3 xdr jf; tfjdr - ^V 2 - n)if> a - V44^ 



The indices a = 1, 2 refer to the two components f, I, of spin. The Grassmann- 
valued Fermi fields if), ifr , are to be taken antiperiodic under r — > r + (3. 

A positive value for V corresponds to an attractive interaction between 
the particles. Given an attractive interaction, and a low enough temperature, 
the system should be unstable with respect to the onset of superconductivity. 
To detect this instability we introduce the ancillary complex scalar field A 
which will become the superconducting gap-function. We use it to decouple 
the interaction by writing 



Z = [ d[ifj]d[4]d[A]d[A*] exp - f d 3 xdr { ]T ipl(d T - —V 2 - 

J JO ta=l ^ m 



-A*V> 2 V>i-AVM + ^|A| 2 }. (2) 

The equation of motion for A shows us that A = Vip2' l l ) i- 

Taking note of the anticommutativity of the Grassmann fields, the quadratic 
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form in the exponent can be arranged matrix 





We may now integrate out the fermions, obtaining the functional determinant 
of the Bogoluibov-de Gennes (BdG) operator 

The BdG operator acts on two-component Nambu spinors obeying ^(x, r + 
/3) = -V(x,t). 

Provided we can ignore quantum and thermal fluctuations of the gap 
function A(x) = |A|e* e , we can perform the integral over A by simply re- 
placing it by its saddle-point value. In this way the euclidean time effective 
action for a BCS superconductor is given as 

r(A) = -lnDetS + yci 3 a;cit^|A| 2 . (5) 

In this work we will be concerned with the situation where the gap-function 
A is time- independent. In this case the functional T reduces to (3 J 7 (A) where 
T is the free energy. 

It is not possible to evaluate T exactly for an arbitrary gap-function, but 
when A is much smaller than E F and varies slowly on the scale of the Fermi 
wavelength, then only momenta near the Fermi surface are important. In 
this case it is reasonable to replace the single-particle energies of the normal 
metal by the linearized form 

e(k)=v F (\k\-k F ), (6) 

and to approximate k 2 d\k\ by k 2 F d\k\. Having done this, then with no further 
approximation, we may then use the results from @ to write 

r(A) = -2nv F N(0) /^c/ 2 x ± lnDet(9 r -w F a 3 (n-V) + |A|a ie - i,J3e ). (7) 
J Aix 

What is happening in equation ([7|) requires some explanation. Firstly 
N(0) is the density of states at the Fermi surface 

. . mk F , . 

mo) = -^f , (8) 
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the symbol n denotes a unit vector, and Vp = kp/m is the Fermi velocity. 
Having chosen a point k = nkp on the Fermi surface, we decompose the 
coordinate vector x into a part parallel to n, and a part perpendicular 



X = X\\U + Xj_. 



(9) 



Then, for fixed xj_, we compute In Det (d T +H) where 7i is the one-dimensional 





Andreev hamiltonian 



H = —ia 3 v F (n ■ V) + |A|<7ie 



-2(73(7 



-iv F (n-V) A 
A* +iv F (n-V) 



(10) 



This Dirac-like operator appears as an approximation to the full hamiltonian 
when we restrict our attention to Nambu spinors of the form x( x ) ex P ikpn ■ x, 
where xi x ) is slowly varying. Here the spinor \ and the gap-function A are 
to be regarded as functions of xn. Notice that the derivative n • V coincides 
with d x ,. . Having evaluated In Det (d T + 7i) , we then integrate the result over 
the family of parallel lines parameterized by xj_. Finally the / performs 
an average over the directions n. 

If we are given an expression for In Det (d T + 7i) in the form 



In Det (d T + H) 



F(x\\, x.±)dx\\dt, 



(11) 



then the integrations over xn and xj_ combine to give an integral over ev- 
ery point in three-dimensional space. To illustrate this, consider a time- 
independent A at zero temperature. To second order accuracy in the gradi- 
ent expansion the part of the free energy coming from gradients of the phase 



of the order parameter is [16 



[^((n-V)6) 2 dx {l dt. 



Using 



dQ n 

4vr 1 ] 



and remembering that kf — mvp, the effective action becomes 

1 



r(A) 



K f 

3tt 2 J 8 



5m 



■(vey<rxdt. 



(12) 



(13) 



(14) 
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This is not unreasonable. Recognising that 

k 3 

& = (15) 

is the number-density of the electrons, and that the superfluid velocity is 

v, = iv«, (16) 

we can rewrite ([14]) as 

r(A) = / ^p s mv 2 s d 3 xdt . (17) 

Since everything is time- independent the t integral is trivial, and from T(A) = 
PJ-'(A) we have 

F(A) = J ^p s mv 2 s d 3 x. (18) 

Thus this part of the free energy coincides with the kinetic energy of the 
superflow. 



3 The Resolvent Equation 

When 7i is time-independent we have 

oo 

In Det (d T +H)= Y, ln Det ( iuJ n + H), (19) 

n=— oo 

where u n = 2ir(n + the fermionic Matsubara frequences, are the eigen- 
values of d T . 

Our task therefore is to obtain the the Fredholm determinant Det (TC—E) 
involving the Andreev hamiltonian operator 

U = -ia 3 d x + |A|a ie - iCT3e , (20) 

and where E is a general complex number. (We have set Vp = 1.) Such 
determinants are usually evaluated by use of the variational formula 

SlnDet (H-E) = Ti{{7i- E)- l 5{7i- E)}, (21) 
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where {7i — E) 1 denotes the Green function, or resolvent, 

G aP (x, y, E) = (a, x\(H - E)~ l \(3, y). (22) 

The Green function can be regarded as an infinite dimensional matrix in the 
variables x, y, and as a 2 x 2 matrix in the Nambu spinor space on which 
the <7j matrices act. The notation TrA includes an integration over the x,y 
labels as well as a conventional trace (to be denoted by tr ) over the spinor 
indices. 

Provided E is not an eigenvalue of 7i, the Green function exists and obeys 

{-ia z d x + \A\a ie - 1 ^ 9 - E) a ^{x, y, E) = 5{x - y)5 a p. (23) 

If V'a and ipa are solutions to the equation 

(-ia 3 d x + \A\a ie -^ e - E)i/> = 0, (24) 

satisfying suitable boundary conditions to the left and right respectively, then 

G a p(x,y,E) = ^i)%{x)i)%,{y){<Ti)pp for x > y 

= wV^(z)V#(2/)<>iV/3 for y>x. (25) 

Here 

W = W^ L , = -t(^ L ) T a 2 ^ R = ~ (26) 

The expression W(ip L , ifj R ) is independent of x for any tf) L , i/j r that are both 
solutions of (TC-E)ip = 0. W plays a role analogous to that of the Wronskian 
in the theory of scalar linear differential equations. In particular W vanishes 
identically if and only if ip L and ip R are linearly dependent. Such linear 
dependence signals that E is an eigenvalue. 

When making the variation 5{7i—E) in (M) we will only consider changes 
in E and A, so we will require G(x, y, E) only at the points x = y. Now 
strictly speaking, G(x,x;E) is not well defined, since the jump condition 
implied by ( p3~D is 

G(x, x + e) — G(x, x — e) = ia 3 . 
Fortunately this ambiguity does not affect 51nDet (TC — E). This is because 

tr[a 3 5( A f A ^)]=tr5(jf, ^)=0. 
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We may therefore safely define the symbol G(x, x, E) to mean the average of 
its two limits [QEll 



G{x,x) = f lim -(G(x + e,x) + G(x,x + e)). (27) 

e-^0+ 2 

With this understanding consider the matrix- valued function g(x) = 
2iG(x,x, E)a 3 . In components this is 

g(x) af3 = ± (v^< + faW (28) 

From the equations obeyed by ip R an d 4> L we immediately see that g obeys 
the equation 

d x g=[M,g], (29) 

where 

M = -ia 3 (\A\a ie - ia:ie - E). (30) 

The reason for the appearance of 03 in the definition of g is that cr 3 appears 
in the coefficient of d x in the Andreev operator (|20|). If we pre-multiply the 
Andreev operator by — io~3, it becomes — d x + M, where the coefficient of 
d x is now proportional to the identity matrix. The Green function for this 
modified operator is G(x,y, E)a 3 , and g(x) is, up to the factor of 2i, the 
diagonal element of this Green function in x space (it is still a matrix in 
Nambu spinor space). 

The equation (|29|) can also be written as a commutator 

[-d x +M,g] = 0, (31) 

where we regard g(x) as an operator acting on the space of vector valued 
functions by pointwise multiplication. 

In the superconductivity literature ( |2"9"D is usually called the Eilenberger 
equation and g the Eilenberger function ||13|| . It lies at the heart of the quasi- 
classical theory of superconductors [|nj. Equation ( |2"9"D is also the analogue, 



for matrix first-order differential operators, of the Gelfand-Dikii equation 



[15|. Like the Gelfand-Dikii equation, it has extensive applications in the 



theory of integrable dynamical systems Because we wish to stress that 
( |29|) is an exact property of the diagonal element of the resolvent of any 
first-order matrix differential operator of the form — d x + M (rather than an 
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approximate property of the partially integrated Gor'kov Green function) we 
will often refer to it as the resolvent equation. 

Using the definition of W, it is clear that the particular solution of the 
resolvent equation given by equation fl2"8|) satisfies the condition tig = 0. 
Similar algebra shows that 

0V = 1> L , 9 2 ^ R = *P R - (32) 

Since the two solutions are linearly independent and span the space of Nambu 
spinors at each x, this implies g 2 — 1. This condition is sufficient to determine 
the solution^ 

Given g as a function of E, it is possible to recover lnDet (TC — E). We 
simply observe that 

^ In Det (H-E) = -Jdx ^tr (ga 3 ) } (33) 

and integrate with respect to E. 

To obtain a gradient expansion of g we introduce a parameter z (which 
will ultimately be set equal to 1) and rewrite (|29| ) as 

d x g = [zM,g]. (34) 

If we seek solutions in the form 

00 

g = Y J 9nZ~\ (35) 

n=0 

we obtain the deceptively simple recurrence equation 

d x g n -i = [M,g n ). (36) 

At first sight this determines g n in terms of a derivative of g n -i, so the n-th 
term contains n gradients of the gap parameter. 

lr The normalization g 2 = 1 is specific to the resolvent of the Andreev operator on an 
infinite line, or on a finite interval with self-adjoint boundary conditions at each end. When 
periodic boundary conditions are imposed on an interval of length L, the normalization 
becomes g 2 = — 1 cot 2 kL/2, where k is the Bloch momentum corresponding to E in the 
periodically extended problem. 
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Unfortunately (^) is not quite as innocent as it looks. The n-th recur- 
rence relation only immediately determines g n up to the addition of terms 
that commute with M. All is not lost however. Because M is a traceless 
2x2 matrix, the only undetermined term is in fact proportional to M - 
and it is possible to find this term from the next equation in the series 

d x g n =[M,g n+1 }. (37) 

Taking a trace shows that tr Md x g n = 0, and this contains the information 
we need. Using this method we have evaluated the gradient expansion up 
to eighth order in the gradients fll5l . Because extracting the term cx M 



requires an explicit parameterization of the space of 2 x 2 matrices, the final 
expressions are given in terms of gradients of |A| and V# separately. It is 
not at all obvious how to assemble these very lengthy expressions into more 
compact polynomials in derivatives of A itself. In the following sections we 
will descibe an alternative approach that does not require dissecting g into 
parts parallel and perpendicular to M, and indeed does not require us to find 
g directly at all. 



4 A Useful Identity 



In what follows we will concentrate on evaluating the determinant in the case 
where TC is defined on the infinite real line R. For convenience we will assume 
that all x dependence of the gap function A(x) takes place in a finite region 
Qq, and that to the left and right of Qq lie regions Ql and respectively, 
in which A takes constant values A^ and Ar. 

What we do next is motivated by the discussion of first-order matrix 



differential operators in |I2 |. We observe that because g is a 2 x 2 traceless 
matrix, it has distinct eigenvalues and is therefore diagonalizable. Because 
g also obeys the equation g 2 = 1, we see that these eigenvalues are ±1. 
There therefore exists a 2 x 2 matrix <p{x) such that g(x) can be written 
g(x) = (j){x)a ? ,(j)- 1 (x). 

As mentioned earlier the resolvent equation can be writtenP] 

[-d x +M,g] = 0. (38) 



2 We hope that it will be clear from the context when a symbol such as d x g refers to 
the derivative g' , and when d x g is an operator acting on functions according to (d x g)ip — 
g'(p + gcp'. 
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Substituting g = 4>a$4> 1 into fl3"B|) gives 

0[0- 1 (-^ + M)0, O r 3 ]0- 1 = 0, (39) 

and this implies that 

[<j ) - 1 (-d x + M)<f ) ,a 3 } = 0. (40) 

Consequently, if we define a matrix A by 

r 1 (-d x + M)<p = -d x + A, (41) 

then A = _1 Af0 — commutes with cr 3 and so must be diagonal. The 

matrix <p(x) is sometimes called a dressing function. 

We now derive a useful identity that connects A with the determinant we 
are wish to compute. We claim that 

tr{g8M) = 5ti (Aa 3 ) - ^d x tr {(jx?^- 1 - S^- 1 ) 

= 5L + d x a. (42) 

(Here and elsewhere 5<p~ 1 denotes 5(<j)~ l ) and not (<50) _1 .) 

We call the two quantities L and a because of the analogy with classical 
mechanics. In mechanics we have 

(^-mqi - ij^J Sqi = 5j2 (^ m i<ii - V^j + dt (- m ^Mi\ ■ ( 43 ) 

Here the first term on the right-hand side is the variation of lagrangian, while 
the term in the total derivative generates the symplectic form which deter- 
mines the hamiltonian structure. The left-hand side contains the functional 
derivative of the action, and gives the equations of motion via D'Alembert's 
principle. Dikii || exploits this analogy when discussing the stationary equa- 
tions of the integrable hierarchies. 

Notice that the second term on the right hand side of (|42|) can equally 
well be written as — c^tr {(j)a 3 5{(f)^ v )) or as d x ii (S^a^cj)^ 1 ). This is easily seen 
by using (j>4>~ 1 = 1. 

To establish ( f4"2Tj we rewrite ([IT]) as 

-d x (j) + M(f) = 0A (44) 
d,^- 1 + cj)- 1 M = Acj)- 1 (45) 
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From these we have 

A = f -djcj) + SM<j> + M5<f> - 50A - (f)5A = (46) 

B = f djcj)- 1 + 4r x 5M + 5c))- 1 M - A50" 1 - 5A0" 1 = (47) 

Thus 

= ti (Aas^ 1 + (f)(j 3 B). (48) 

Expanding this out, using the cyclicity of the trace, the equations _1 M = 
A0 _1 — d x (j)~ l and M0 = 0A + <9 X 0, together with g = 0cr 3 0~ 1 , we find (|42|). 

To see the utility of ( f42|) observe that from g = 2iGa 3 and <5M = 
— zo"3<5(7i — i?) we have 

/oo 
dxtr(g(x)8M(x)) = 2Tr (G8H) = 28 In Det (W - E) . (49) 
-oo 

Thus, if it is legitimate to ignore the boundary terms coming from the total 
derivative, we have 

roo \ 

In Det (H-E) = / -tr (Aa 3 )cix + constant. (50) 

J— oo 2 

Here the "constant" is a quantity unaffected by local changes in A(x). It 
will, however, depend on the asymptotic gap functions A LjR . 

In general the boundary terms cannot be ignored, but it is possible to 
choose so that they are zero. We have this freedom because requiring 
that g(x) = 0(x)ct30 _1 (x) does not uniquely determine 0. Indeed the re- 
placement — > <px with xi x ) an y invertable diagonal matrix leaves g(x) 
unchanged. Such a "gauge transformation" does however alter A. Under the 
above substitution we have 

A^A-rtx. (51) 

By making such a gauge transformation we can transfer any of the contri- 
butions to tr (g8M) from A to the total derivative term, or vice versa. Two 
extreme choices come immediately to mind. 

1. Select so that A = 0. In this case all the contributions to the deter- 
minant come from the boundary term. 
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2. Select so that the boundary terms constribute nothing to the deter- 
minant. 

The latter choice is possible because the Green function G(x, y, E) is in- 
sensitive to the values of A at large distance from x and y. We have assumed 
that there are asymptotic regions £Il,r where A(x) becomes independent of 
x. Once x is well inside either of these regions then g(x) will settle down 
to a constant value depending only on the asymptotic values ^l,r- We can 
therefore select a <fi depending only on A^ to diagonalize this constant ma- 
trix. Once we have done this then 0(±oo) will be unaltered by variations of 
A(x) taking place in Qq- The determinant is then obtained entirely from A. 
This is the route we will take to compute the gradient expansion. 



5 Connection with "Shooting Method" 

Before we derive the gradient expansion, we should point out that our first 
choice of 0, the one leading to A = 0, is equivalent to the so-called "shooting 
method" for evaluating the determinant |Tj|. This asserts that, in a suitably 
interpreted sense, the determinant is proportional to the inverse of the trans- 
mission coefficient for scattering off the spatially varying A(x). To see this, 
observe that if we set A = 0, so that 



<f ) - 1 (-d x + M)<f ) = -d x + 0, 
then we are requiring the columns of to be solutions of 

(-d x + M)^ = 0. 

Equivalent ly 



A candidate for 



(-ia 3 d x - 
is therefore 



A\a ie - ia3e - E)ip = 0. 



(52) 

(53) 
(54) 

(55) 



where the colums are composed of the ip L and ip R solutions used in construct- 
ing the Green function. The inverse is 



1 

W 



^2 -^1 



R 



(56) 
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where, as before, W = W(ip L , ip R ) = ip R ip2 ~ V'lVi'- 

A short computation shows that 0<730 _1 = g(x), so we have guessed 
correctly. 

Now we write down ip L and ip R in the asymptotic regions, which is the 
only place where we will need them. Let us assume that E is real and that 
E 2 > |A^| 2 , |A^| 2 , so we have scattering solutions. 

When A = |A|e* e is constant, a plane-wave solution to 

(-ia 3 d x + \A\a ie - ia3d )^ = Eip (57) 

is given by 

ijj = u(k, E, A)e ikx = 

Here E 2 = k 2 + |A| 2 . 

We now introduce transmission and reflection coefficients, t^L and r R L 
and define 

^ R ( x ) = u (k, E, A L )e ikx + r L u(-k, E, A L )e~ ikx x e tt L 

= t L u(k,E,A R )e ik ' x xeVl R (59) 



(58) 



il; L ( x ) = u (-k,E,A R )e' ik ' x + r R u(-k',E,A L )e ik ' x x E tt R 

= t R u(k, E, A L )e~ ikx xen L . (60) 

The apparently perverse appearance of the subscript L on the reflection and 
transmission coeficients in ip R (and mutatis mutandis in ip L ) is supposed to 
indictate that the incoming wave was incident from the left (right). The 
wavenumbers k and k' are not necessarily equal since we do not assume that 
|A L | = |A fl |. 

Substituting these expressions into W(ip L , ip R ) we find that 

W(^j l (x), iP R (x)) = -i[^ L fa^ R = -2kA L t R xett L 

= -2k'A R t L x E Vl R . (61) 

Since W is x independent, we must have 

kA L t R = k'A R t L . (62) 
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This reduces to ti, = tn in the particular case that Ax, = Ar. 

In order to use if) and ip L in the Green function we require that k, k', 
have a positive imaginary part so that if>R tends to zero as x — > +00, and 
similarly if>i tends to zero as x — > — 00. If we assume that both the real part 
of E and the real part of k are positive (so that if) and ip L correspond to 
a real scattering process, and G to outgoing waves) then this requires the 
addition of positive imaginary part to E. If we wish to evaluate G, and from 
this Det (Ti, — E), for E below the positive real axis then we must replace k 
by —A; in the above wavefunctions, so that the resulting negative imaginary 
part of k still leads to convergence. As expected, this means that both G and 
Det (Ti — E) have a branch cut discontinuity across the real E axis whenever 
asymptotic plane-wave solutions are possible. 

Now we use these functions to evaluate tr (50cr 3 _1 ). Near x = —00 we 
have if) R — > u(k, E, Ai)e tkx which is large, while ip L = t^u(k, E, Ax)e~ lfca: is 
small. The only expressions occuring in tr (<50a"30 _1 ) involve their product 
which is finite. Varying A in Qq changes tn, while leaving u(k, E, Ax) etc. 
unaltered. We find that 



Similarly 



tr (50o- 3 X ] 



St 



L 



(64) 
ti 



Inserting these results into (^) leads to 

In Det (Ti - E) = - In t + constant. (65) 

It does not matter whether we use tx or t^ in this formula because the 
logarithms differ by terms involving E and Axj,x, and these can be included in 
the constant. The constants cancel when we consider ratios of determinants 
of operators with the same asymptotic A's. In other words, if Ti and Ti^ 
are two hamiltonians with the same asymptotic behaviour, then 

Det (Ti — E) tf(E) tf '(E) 



Det - E) t L (E) ' tjt(E) 



(66) 



Notice that the (analytic continuation of) t and r become infinite as E 
approaches the energy of a bound state. The wavefunctions ip R and ip L be- 
come proportional to one another, and decay exponentially as \x\ — > 00. The 



16 



resultant vanishing of the determinant is consistent with the interpretation 
of the ratio of the Fredholm determinants as a regularized version of 



(K-E)_ 



where the E n are the eigenvalues of H, and the the eigenvalues of H^. 
For E in the continuous spectrum the zeros and poles coming from the eigen- 
values of H and merge to form the branch-cut noted above. 

6 The Gradient Expansion 

To find the recurrence relation for we multiply the matrix M by z as before. 
At the end of the calculation, we may set z — 1. Thus again we are looking 
for a solution of 

[-d + zM,g] = (68) 

as a power series in 1/z. To find the free energy expansion, we do not need 
g itself but only and A that satisfy 

zM(f) = <p' + <f)zA, (69) 

and are given as power series in 1/z. 

For the homogeneous superconductor, only the zeroth order term in A 
will be non-zero. We find it from the zeroth order expression for g. Since g 
is traceless, squares to the identity, and commutes with M, we must have 

9o = jM, (70) 

where ( is defined by M 2 = ( 2 1. Since 

M = (?r. ~ iA ) (7i) 



(we have now replaced — E by the Matsubara frequencies) we have 

C = V^ + l^)! 2 - (72) 
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From (0) and fl70|) we find 

5tr(a 3 A )=tr(g 5M)= 25(, (73) 

so 

tr(A a 3 ) = 2C, (74) 

and 

Ao = (75) 
Hence, we look for A in the form 

A = C°"3 + — + -| + • • ■ (76) 

Comparing coefficients of in (|6"5|), we obtain the equation for the 0-th order 
in 0: 

M0 O = o Ca 3 (77) 

To obtain a recurrence relation for higher-order terms in 0, it turns out that 
0o has to be factored out from the expansion, i.e., we look for in the form 

= 0o(l + ^ + $ + ...). (78) 

Now we can go on to calculate the first-order terms in ( |7B] ) and ([7^) by 
comparing the coefficients in fl69"|): 

M0o0i = 0o + <Mi + 0o0iC^3 (79) 
Multiplying by 0q 1 on the left and using ([77|), we obtain 

C[a 3 ,0 1 ]-A 1 = o ' 1 0o- (80) 



In general, comparing the coefficients of z n in (|69"D, we get the n— th order 
relation 

C[d 3 , n+ l] - A n+1 = 0^ + 0o V O 0n + 0nAl + 0„-lA 2 + . . . + 0iA n . (81) 

The in + l)th terms in the expansions of and A are thus given by all 
the previous terms. Moreover the first term on the left-hand side of (|8l|) is 
perpendicular to a 3 , whereas A contains a multiple of <r 3 and a multiple of 
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1, since it commutes with 03. In other words, the first term on the left-hand 
side contains o\ and 02, whereas the second one contains 03 and 1. Hence, 
we obtain the (n + l)-th terms in the expansion of <fi and A by calculating 
the right-hand side of (^T|) (in which everything is known from the preceding 
steps in the recursion), and splitting it up into the transverse and longitudinal 
parts relative to 03. 

There is a minor hitch: (^l|) does not determine the part of </> n +i that 
commutes with 03. Adding a term to n +i that commutes with 03 does 
change the right-hand side of all subsequent recurrence relations. Such a 
term may, however, always by removed by changing our choice of the <fi that 
diagonalises g. This will change the free energy density by a total derivative 
of terms local in the gap. It will not affect the total free energy. 

We may therefore proceed to find and A. Unfortunately this will lead 
to the same lengthy expressions we found in [jnj by evaluating the recurrence 
relation for g. The problem is that there is no simple way to write the matrix 
4>q Vo' in terms of M itself. We will therefore have to choose an explicit basis 
for the matrices and work with components. 

After a lot of thought we realized that is far more convenient to replace 
the arbitrarily chosen matrix 03 in the decomposition g = 0030~ x by the 
local matrix M. To do this we transform 

-»• $ = 00 o = 1 + + + . . . 

Z Z l 

= i + — + -£ + ... 

Z Z A 

A^A ee = M+ ° Al00 " 1 + 0oA2 /°" 1 + ■ ■ ■ 

z z 2 

M + — + § + ••• ( 82 ) 



z z 

Then L = tr (A03) = tr (AM/(). By using 



00^00 1 = K ~ 0000 1$ n + $n0' O 0O \ 

we can rewrite fl8lD as 

[M, $ n+1 ] - A n+1 = $' n + $„0' Q 0o 1 + ^nAi + . . . + $ n Ai. (83) 

The advantage of this expression is that the two factors of $ n in the 
second and third terms on the right hand side are both on the left. This 
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enables us to get rid of the awkward expression 4>' 4> 1 by factoring out $ n 
and using the recurrence relation for n = 



We obtain 

[M, $ n+1 ] - A n+1 = & n + $„[M, $ x ] + $ n _!A 2 + . . . + $ n Ai. (84) 
We need, however, to start the recurrence by obtaining $ x in some way. For 



this we can use the Eilenberger equation 
the expansion of g in 1/z 



) . From 



?o"3<; 



C 



Af 1 



M 



and (|77), we get 



+ . . . 



. 9i 
9o + — + --- 

z 



The coefficient of z in (pq) gives 



that is, 



M, 



* M 



The double commutator of matrices A, B, C can be written as 
[A, [B,C\]={{A,B},C}-{B,{A,C}}. 



(85) 



Here the braces denote an anti- commutator . For three traceless 2x2 matrices 
this reduces to a form of the vector triple-product identity 

[A, [B, C}} = 2tr (AB)C - 2tr (AC)B. 

We now choose $i with no longitudinal components, i.e. tr $i = tr ($jM) : 
0. Using trM 2 = 2( 2 , we get 



$ 1 = - 



1 1 (M 



(86) 
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We see the natural appearance of Mj ( as well as of the derivative multiplied 
by Hence, it is natural to divide fl34|) by (, and to define 



N = 



M 

T 



—r— A Ar 1 Ai 1 A 2 

»t ^ A 2 

2; 2r 



and 



D = la. (87) 

The symbol .D still behaves as a derivative in the sense that it is linear, and 
obeys the Leibnitz rule 

D(AB) = (DA)B + ADB. 
In terms of these new variables, we can rewrite ( [34]) as 

[N, $ n+1 ] - A n+1 = D1> n + $ n [A r , $j] + $ n _x A 2 + . . . + $! A n , (88) 

and 

L = Ctr(AJV). 

The formula fl5B| ) leaves undetermined the part of $ n +i that commutes 
with AT (just as fl8lD did not determine the part of 4> n +i that commuted with 
(T3). Again it is convenient to choose the $'s purely transverse, i.e., 

tr $ fc = tr ($ k N) = 0, for all k > 1. 

With this choice we get a recurrence relation for the $'s only, and a very 
simple expression of L in terms of the $'s. 

To see how this comes about, let us think of N as being 03 (this can always 
be achieved at a given point by a global rotation which does not change L), 
then the $'s contain o\ and 02, whereas the A's contain 1 and 03. Hence, 
among the terms in d88|) , A n+1 and $ n [AT, $1] contain 1 and o" 3 only, [A/ - , $ n +i] 
and all the products of the type §jA k contain o\ and 02 only. Finally -D$„ 
potentially contains all four matrices. We can separate the transverse and 
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the longitudinal part by taking the commutator and anticommutator of N 
with both sides of (|88|). From the commutator we get 

[N, [N, $ n+1 ]] = [N, D<t> n ] + [N, ^ n _{K 2 + ... + $ x A n ]. (89) 

Using the double commutator formula (|85[) and its special case for traceless 
matrices, this leads to 

[N, [N, $ n+1 ]] = 2tr (7V 2 )$ n+1 - 2tr (N<t> n+1 )N = 4$ n+1 , (90) 

where the second term in the middle expression is zero due to transversality. 
For the last term on the right-hand side of ( |8"9"D we use the formula 

[A, BC] = {A, B}C - B{A, C} 

to get 

[N, $ fe 4-] = {N, $ fc }4. - $ fe {iV,A,}- (91) 
Since {a 3 , 01,2} = 0, the first term on the right-hand side of (|9~ID is zero. So 



we can rewrite ( |89| ) as 

4$ n+1 = [N, D$> n ] - ^{N, A 2 } - ... - $i{iV, AJ. (92) 

We see again the usefulness of introducing N,A, and D. 

Notice the appearance of the anticommutators {N, A,} on the right-hand 
side of (0). Due to transversality of the $'s, the A's themselves are elimi- 
nated. We can relate these anticommutators to the $'s by taking the anti- 
commutator of iV with both sides of Q58]). 

-{iV,A n+1 } = {N, D$ n } + {iV,$ n [iV,$ 1 ]} = 

= D ({N, $„}) - {DN, $„} - i {N, $ n [N, DN}} . (93) 



The first term on the right-hand side of ( |9~3"D is zero just as in fl9T|). In the 
last term, we use 

= Dl = D(N 2 ) = N ■ DN + (DN)N 

to see that 

[N, DN] = 2N ■ DN. 
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Then we use the formula 

{A, BC} = {A, B}C - B[A, C] 

to get 

- ^{N, $ n 2iV ■ DN} = ~{N, ® n }N ■ DN + U n [N, N ■ DN}. (94) 

The first term on the right-hand side of ([)4]) is again zero; for the second 
one, we need to evaluate 

[N, NDN] = N 2 DN - N(DN)N = DN + N ■ N ■ DN = 2DN, 

where we used again the anticommutation of N and DN. So we can rewrite 
(H) as 

-{iV,A n+1 } = -{DN,$ n } + $ n -DN = 

= -(DN)-$ n . (95) 

We introduce the symbol 

C k = {N,A k } 7 (96) 



and substitute from (^) to (^) to get a recurrence relation for the $'s only: 
4$ n+1 = [N,D$ n ] - ^(DN)^ - ... - ^(DN)<t> ni . (97) 
We can then simply calculate the terms in the expansion of £ = {N, A} as 

£ n = {DN) ■ (98) 

This gives us the desired lagrangian as 

C 

L = |tr£. (99) 

The last three equations provide all the machinary needed to compute 
the first few terms in the expansion of L. We start with 

$! = --DN 
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from 



). Then 



4$- 



[N,D^] = 
-~[N,D 2 N], 



so 



$0 



16 



[N,D 2 N]. 



In the third order, we have 



4$, 



[N,D$ 2 ] ~ $i(.DjV)$i 



16 
1 



1 



N,D[N,D 2 N}\ - ^ DN ^ 
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1 r 
16 



1 



N,[N,D 3 N] ~-(DNY 
J 16 



N, [DN, D 2 N] 

We now use the formula (^) for the double commutator and get 

$ 3 = - ^{{n,dn},d 2 n} + ^{dn,{n,d 2 n}}- 

- h {{*> N h D 3 N} + 1 {N, {N, D 3 N}} - ±-{DN) 



These expresions can be simplified. We begin by evaluating the anticom- 
mutators 



2N 2 = 21 


D ({N, DN}) - {DN, DN} = -2(DN) 2 



{N,N} 
{N, DN} 
{N, D 2 N} 
{N, D 3 N} 



All the anticommutators are proportional to the unit matrix (they are squares 
of traceless 2x2 matrices), so 



D [{N, D N}j - {DN, D N} = -2D ((DN) ) - D {(DN) 
-3D ((DN)' 



±{DN, (DN) 2 } - D 3 N} ~^{N,D ((DN) 2 ) } - ^(DNf 

1 3 
(DN) 3 D 3 N -ND ((DN) 2 ) . 



32 
5 
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To obtain the gradient expansion of the free energy up to the fourth order, 
we need to calculate C 2 and £4 from (|98|): 

C 2 = (DN)& 1 = — ^(DN) 2 

A = (DN)<S> 3 = -^(DN)*-±(DN)-(D*N)-^(DN)-N-D((DNf), 
so 

L 2 = -|tr ((DN) 2 



' Ctr ((DN) 4 ) (tr ((DN) ■ (D 3 N)) - —(tr ((DN) -N-D ((DN)' 



128" v v ' I 32" v ' ^ ') 64 
We can rewrite L 4 in a more symmetric way by writing 

Ctr ((DN) ■ (D 3 N)) = CD (tr ((DN) ■ (D 2 N))) - (tr ((D 2 N) 



The first term is a total derivative ((D = d by definition (|87|) ), and is, 
therefore thrown out. The last term in L4 is equal to zero too, since 

tr ((DN) ■ N) = ^tr ({DN,N}) = 0, 

and D ((DN) 2 ) is proportional to the unity matrix. Thus, up to the fourth 
order in derivatives, including even-order terms only, 

L L + L 2 + L 4 = Ctr {l - X -(DN) 2 - ^(DN) 4 + ^(D 2 N) 2 } . (100) 



7 The Free Energy 

The equation (|100|) at the end of the previous section is the principal result 
of this paper. It contains, in an extremely compact form, all the information 
needed to obtain the free energy. Using the results of section 2 this may be 
written 

T = f rf3 x l H - H «! 2 + f d 3 x \^l + 27T N f^f d 2 X ± ^ ld \ (101) 
J 8lY J V J 47T J 
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where is given in terms of L as 

r 

2 



= -Zj2 f dx\\L. (102) 



Since <i 2 x_|_<ia;|| = <i 3 x, we can write 

r f |H- H a | 2 |A| 2 „ dfi n ] . . 

Notice that terms with an odd number of derivatives, and hence an odd 
number of vectors n, drop out when we average over them. This explains 
why we calculated only even terms in the expansion of L. 

We have included a magnetic field in ( |101| ). In our earlier calculations, 
we did not mention the coupling of the electrons to the magnetic potential 
A. However the component A\\ can always be gauged away along the line 
xi = const, by absorbing it into 9. There is therefore no loss of generality in 
our formulae. To insert A we merely replace our derivative of the complex 
order parameter by a covariant derivative, and our derivative of the (real) 
magnitude of the order parameter by a plain gradient. 

We would like to compare our calculation with that of reference Q as this 
seems to be the only place where the fourth-order terms have been written 
down explicitly. 

The expression given for the free energy in || is, in the notation of that 
paper, 



T = J dh[ |H + N (3- 2 [P 2 w + ^) 2 {£|Ox| 2 + ^(V|x| 2 ) 2 } + 

+ l(^/?W|o 2 xl 2 + ^[||o x | 4 - |Ox| 2 (V 2 |x| 2 ) + ^(V 2 |x| 2 ) 2 ] + 

+ /'[J(V 2 ix| 2 )(V|xi 2 ) 2 - i|Ox| 2 (V| X | 2 ) 2 ]}]}. (104) 

Here x = is the dimensionless order parameter, and O = V — 2ieA 
is the covariant derivative. The homogeneous part of the free energy comes 
from u>(|x|), while g is a function of \x\ 2 that we will identify later (g is not 
to be confused with the Eilenberger function). The primes on g denote the 
derivative with respect to \x\ 2 , i-e. 

' = dg 
9 ~ d\ X \ r 
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The directional averaging / is not stated explicitly in this formula, but is 
to be understood. In other words the product of two expresions AB where 
A, B are either of the two vector operators O or V is defined as 

AB = ^SijAiBj, (105) 
and the product of four factors as 

ABOD=i(^ +Wj+WW H (106) 

We have inserted a question mark over the equals sign in ( |104j ), because we 
believe that there is an error in the fourth-order terms. 

We will now expand our expression so as to write it in Tewordt's form. 
During the calculation, we will keep the primes as derivatives, and only at the 
end we will replace them by the covariant derivatives or gradients according 
to the rule 

x' - Ox 

(M 2 )' - V|x| 2 - (107) 

Note that the prime over the order parameter means dj dx, whereas the prime 
over g means d/d\x\ 2 - 

We introduce the dimensionless quantities 

M=[3M=(( 2m+ y 71 ) (108) 

V IX -(2m + l)7i J v ' 

and 

£ = 0C= v/(2m+l)% 2 + | X | 2 , (109) 

so 

N=^-. (110) 

Since Tewordt writes out the Fermi velocity v p explicitly (in our calculation, 
we set Vp = 1), the dimensionless derivative D equals 

D = V -^d= V -^d. (Ill) 
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From ( |103| ), the second-order term in T is 

= N (3-^(v F (3) 2 J d 3 r £ f e^jtr (DN) 2 , (112) 

where we pulled out the prefactors that appear in the Tewordt formula. Using 
( |TTUD and (pip , we get 

1 -DN l(M ^' M ' Mi ' 



v F p eve; e 2 e 3 

To obtain the result in Tewordt 's format, we need to trade the derivatives of 
e for the derivatives of \x\ 2 - From (|109|) , we see that 

e' = ^p-, (us) 

so 

1 M M(\y\ 2 )' 

DN= J -4- J ^ ] . (114) 



Hence, 

1 o M' 2 (\y\ 2 )' M 2 (\y\ 2 )' 2 



From ( |108| ) we see 

7w 2 = e 2 i (H5) 

{M,M'} = (Mf=(\ X \ 2 )'l (116) 

-M' 2 = (H7) 



so 



1 (DN) 2 = I - ^2 I 1. ( lis) 



M?r y v e 4 4e 

We now use fll07Q to replace the derivatives, and obtain 

H = N^-i^fj^l i^f - ™f\ . (119) 
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Comparing the first term to ( |104| ), we see immediately 

We will also need higher derivatives of g with respect to \x\ 2 - Using ( |109| ), 

we get 

3 _ 1 
q = 7i > — 

15 „ 1 

3 = i^Sr 

9"' = -^Es ("I) 



Hence, 



in agreement with ( pL04j ) . 

The fourth-order term is equal to 



^2 = N (3^(v F (3) 2 J d 3 r (g\0 X \ 2 + {(V| X | 2 )j (122) 



We can calculate the first term immediately from ( |118| ) 

1 (DNr= ((^.xYW)\(jgp, 1 (124) 



To evaluate the second term in ( |123| ) , we need to go back to ( [L14j) : 



1 _ D 2 N = UM' M(\ X \ 2 )'\ 



AT _ 3 ^(lxl a ) ; , M I (\x\ 2 f _ M'" ' 

2 e £ 7 -'c" 
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(we used again ( |113| ) to obtain the second line), so 

a,2 n KA ,2n ,ox/2 



-(D 2 N) 2 



M" z 9 M' 2 (\x\^ 



,, (Ixl 2 )' 2 (M 2 )"' 



M 



+ {M",M} 



(Ixl 



2\> 2 



10 



-{M',M}(\ X \ 2 ) 
2£8 



2£ 5 

(1x1 s 



j2 



-.12 



(\x\ 2 f 

2e io 



From (TO , (1TT5D-(|TT7D, we see that 



//2 



{A4",.M'} 
{.M",.M} 



xV"i 

(-M' 2 )'l 

{A<',A^} , -2A^' 2 = ((| x | 2 ) /, -2 x V / )l ) 



so 



-(D 2 N) 2 



{v F (5) 



x"r" i xWf _ i (ixi 2 ) /4 _ M) 

2^- * ta 4 £10 £8 



Putting the two terms together gives 



iVo/3" 



12 



■/ h U6 32 + 



i6 2 e 



n/1 



ii 



+ 



3 xV" + 3 (|x| 2 ) //2 9 , , M _ 9 (IxrAlxlY _ 3 gWT ' 

r7 + 8 WA J ^ 7 



4 £E 



i6 e 



16 



4 



The two underlined terms do not appear in Tewordt's formula, so we need 
to integrate them by parts (again neglecting the boundary terms): 



i6 2 e 



111 

16 x 8 



(IxP 



f i (lxl 2 ) / 
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91(5 x 8^ ^ \e 



1 HI „ |2V 3/ 1 

1 37 / 2 ,3 lV 37 (\ X \ 2 f(\x\ 2 )" 
316 x 8 1 11X1 j ( 9 j 16x8 £ 9 



and 



glAA ; ^ 7 8 1 £ 7 / g^^ £ 7 16 £ 9 

Therefore, 

^ = No(3 -^J dr ^ 7r {T6-^ + ^2 £5 + 



+ 



,16 £ 7 32 £ 9 4 £ 5 

3 (|X| 2 ) //2 _ 35 (|x| 2 ) /2 (lAf)" _ 15 X 'X*'(\X\ 2 )' 
16 £ 7 16 x 8 £ 9 8 £ 7 



Using ( |121|) and introducing the covariant derivatives and gradients, we fi- 
nally get 

+ /'[4|Ox| 2 (V|x| 2 ) 2 + ^(V 2 |x| 2 )(V| X | 2 ) 2 ]). (125) 

We see that our formula agrees with Tewordt's result except for the pref actor 
in the last term (1/24 instead of 1/4). This is presumably a typographical 
error. We believe our result is correct, because we had previously obtained 
it by two other methods. 



8 Discussion 

While principal result of the paper is the fourth-order term in the free energy, 
as expressed in equations (|100|) and (|125|) , the methods used to obtain this 
term are of interest in their own right. The useful identity linking the dressing 
function <ft with the determinant is one that we have not seen before, and 
the purely algebraic (requiring no integration) generation of the terms is 
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conceptually simpler than other methods of obtaining such series such as 



those we used in ||15|| . In particular the present algorithm produces very 
compact expressions, and is applicable to the computation of determinants 
of any 2x2 first-order matrix operator. 
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